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Inversion techniques are used to infer bottom 
sedimentary characteristics using AN/SQS-53C Reverberation 
Level (RL) data gathered during Exercise LWAD 99-1. 
Analysis was conducted to determine the correlation between 
sediment type and beam RL time series statistics, including 
slope and detrended higher order statistics (standard 
deviation, skew, and kurtosis) . A geographic contour plot 
of the spatial distribution of sediments was generated from 
the deviation of the RL slope of individual beamform data 
from the area-wide average RL slope. The resulting plot 
shows features similar to that of the pre-exercise 
sedimentary characterization map. Higher order statistics 
are found to be non-Gaussian which has signal processing 
implications . 
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I. 



INTRODUCTION 



A. OBJECTIVE 

The objective of this thesis is to develop a real time, 
high-resolution method for determining bottom sediment 
acoustic properties in a shallow water environment. This 
study sought to determine those properties through the 
analysis of the slope characteristics of beam reverberation 
data and higher order statistics (HOS) related to the RL 
beam time series. 

Current fleet anti-submarine warfare (ASW) operations 
are limited in shallow water in part because of their 
dependence on bottom sediment type maps that are based on 
large area averages (low resolution) derived from few 
samples [Ref. 1]. Such maps provide limited or no 
information on the bottom roughness or the sediment 
composition beneath the water-sediment interface. They 
further limit the accuracy of ASW operations because of the 
high spatial variability of the sedimentary composition, 
both laterally and vertically, that is a common feature of 
shallow water continental shelf areas. 

This thesis examines the possibility of using inversion 
techniques to deduce sedimentary properties from beam 
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practical not 
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available for research, but because the RL data is gathered 
by active fleet assets rather than an oceanographic research 
vessel. In this study, the RL data was gathered during the 
Littoral Warfare Advanced Development (LWAD) Exercise 99-1, 
aboard a Ticonderoga class cruiser, using the AN/SQS-53C 
tactical sonar. 

B. SHALLOW WATER ASW PROBLEM 

Since the end of the cold war, anti-submarine warfare 
(ASW) operations have moved geographically from the deep, 
open ocean environment, to the nearshore, shallow water, 
littoral environment. Cold war ASW operations had their 
primary emphasis on locating and tracking Soviet ballistic 
missile submarines (SSBN) on patrol in the open ocean. Over 
the last 10 years, significant geo-political events have 
occurred which have shifted ASW emphasis to the littoral 
environment. With the downfall of the Soviet Union came a 
significant decrease in their SSBN patrols, while 
concurrently new threats arose in such places as the Middle 
East and the Korean peninsula. Smaller countries looked to 
acquire diesel-powered submarines for coastal defense. 
These vessels have become relatively affordable causing 
their numbers and countries of ownership to proliferate. 
The U.S. Navy ASW doctrine has responded to this paradigm 
change by placing increased emphasis on shallow water 
operations [Ref. 2]. 
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Shallow water ASW operations are significantly more 
challenging than those conducted in the deep, open ocean. 
The littoral environment includes ocean fronts and eddies, 
stronger currents, extensive, often steep topographic 
changes, increased sea life, varied bottom sediment 
distributions, and increased shipping traffic. These 
factors combine to create a challenging ASW environment. 

C. REVERBERATION 

One of the most significant factors contributing to the 
complexity of shallow water ASW operations is the spatial 
variability of the sea bottom composition. The reflective 
nature of the overlying sediment exerts a significant effect 
on active sonar reverberation levels (RL) . Reverberation is 
the sum of all random scattering mechanisms in the ocean 
[Ref. 3, p. 237]. In deep water (Figure 1) reverberation 
primarily exists as a combination of two scattering 
mechanisms, volume scattering and sea surface scattering. 
Volume scattering occurs due to scattering from marine life 
within the water column and from inhomogeneities in the 
ocean thermal /density structure. Sea surface scattering is 
due to random reflections (scattering) from a roughened air- 
ocean interface, and is a function of wind speed or wave 
height . 

In shallow water (Figure 2) reverberation from the sea 
bottom typically dominates that due to volume or sea surface 
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reverberation. The contribution of bottom scattering to the 
total reverberation level (RL) is often so severe that a 
target echo that may have been detected in deep water can 
not be detected in shallow water (Figures 1 and 2) . 

D. BOTTOM BACKSCATTERING 

Unlike RL, which is scattering in all azimuthal 
directions, backscattering is only that component of the 
scattered energy that returns to the receiver in monostatic 
(co-located transmitter/receiver ) systems. The intensity of 
the backscattered energy from the sea floor that arrives at 
the sonar receiver is parameterized as the scattering 

strength, 3 =10 log — — [Ref. 3, pp. 238-239]. I scac is the 
b I. 

inc 

intensity of the scattered sound per unit area at a unit 
distance and I iM is the incident plane wave intensity 
impinging on the scattering area. Measurements of RL and 
S b , arising from sea bed backscattered energy, have 

indicated that S b is a function of grazing angle (0) and 

bottom type. Backscattering strength traditionally is 
represented using a Lambert's law relationship where S b = 10 
log (i + 10 log sin 2 (0) . 0 is the grazing angle and fl 

represents the degree of reflectivity of the sediment 
interface. Based on the deep water data of Mackenzie [Ref. 
4] and others, Urick [Ref. 3, pp. 271-273] reports that over 
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moderately rough topography of varying sedimentary makeup 
and small to intermediate grazing angles Mackenzie's 
constant (k = 10 log fi.) averaged about -29 dB, being 

essentially invariant in frequency and grazing angle. More 
recently, shallow water experiments have indicated an 
association of Mackenzie's constant to the bottom sediment 
type and that the grazing angle dependence may be scaled as 
sin(0) over muddy, absorptive bottoms where volume 

reverberation from the upper layers of the sediment may be 
important [Ref. 5]. This finding suggests a direct 
correlation exists between backscattering strength and 
sediment type. With this knowledge, it is possible to use 
inverse techniques to infer the nature of the sediment type 
as a function of its reflectivity from RL data. 

E . INVERSION TECHNIQUES 

The "forward problem" in acoustics are those methods in 
which one predicts the characteristics of the acoustic 
propagation field from a knowledge of the governing oceanic 
parameters (e.g., sound speed structure, sea surface 
roughness, geo-acoustic properties, etc.) through which the 
sound is travelling. "Inverse" (or inversion) techniques 
are the exact opposite such that one infers one or more 
oceanic parameters through the analysis of acoustic 
propagation data. [Ref. 6] 
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Inverse techniques have a wide variety of applications 
including high resolution geo-acoustic modeling [Ref. 7] and 
simulation of broadband signal propagation [Ref. 8] . 

Inversion techniques (IT) are used in this study to 
identify features in the beam reverberation time series from 
which bottom sedimentary properties can be deduced. RL 
features such as the slope of the RL time series and higher 
order statistics (HOS) describing the shape of the RL time 
series are analyzed to identify correlations to sediment 
type. 

F. EXERCISE LWAD 99-1 

The Littoral Warfare Advanced Development (LWAD) 
Exercise 99-1 was conducted from 8 to 15 February 1999 in 
the Gulf of Mexico approximately 200 km WNW of Key West, 
Florida, on the southwest Florida continental shelf (Figure 
3). This test was the ninth in a series of exercises 
performed as part of the LWAD project, a program of shallow 
water exercises designed to identify new technologies which 
have a potential to be transitioned to fleet systems in the 
near future. The LWAD program efforts are concentrated on 
technologies that will aid in the location and 
classification of diesel-electric submarines and mines in a 
shallow water littoral environment. [Ref. 9] 

LWAD 99-1 included several experiments, one of which 
was for the continued development of inversion techniques. 
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The inversion technique experiment had the following 
objectives: (1) record AN/SQS-53C data in a geographically 
variable environment and (2) collect ground truth 
information in order to test monostatic and bistatic 
inversion techniques [Ref. 10]. 

The data examined in this thesis was collected during 
event 31 conducted on 9 February 1999 from 0730Z to 1800Z at 
LWAD 99-1 site 1 (Figure 3). The USS Vicksburg ship's track 
during this event is shown in Figure 4; over 900 pings from 
the AN/SQS-53C tactical • sonar were recorded during the 
event. Event 31 was conducted with a target present; the 
target track history is also shown in Figure 4. 

1. Geo-Acoustics of the South Florida Platform 

The LWAD 99-1 exercise area lies on the southwest 
portion of the continental shelf bordering the west coast of 
Florida. The shelf is composed of a karstic Miocene 
platform that is under Pliocene-Pleistocene and Holocene 
sediments [Ref. 11]. The composition of the surface 
sediment was analyzed by Holmes from a variety of sources 
including sediment core samples and seismic data. Figure 5 
is a depiction of the surface sedimentary distribution at 
LWAD 99-1 site 1 based upon Holmes' analysis and represents 
the most definitive characterization of the upper ocean 
sediment composition prior to the commencement of LWAD 99-1. 
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Howell Hook reef, lying along the western portion of 
site 1 at a depth of 150 m is a natural border between the 
inner shelf and the shelf slope. Howell Hook is composed of 
coral reef material, while immediately to the east is a 
region of inter-reef sediments. The inner shelf, a region 
of recent sedimentary processes (depositional) , contains a 
variety of surficial sediments and is characterized by a 
gentle slope of 0.07°. Farther eastward on the inner shelf 
is a region of foraminif eral sands with sandwaves common. 
Westward of Howell Hook the slope is approximately 0.2° and 
increases through a series of steps until the escarpment is 
reached with slopes as much as 45° [Ref. 12]. 

As part of the environmental characterization for 
Exercise LWAD 99-1, sediment samples were collected. 
However, the results were not available in time for 
inclusion in this study. Preliminary analysis indicates 
that Howell Hook is overlain by a layer of silty sand, while 
the inter-reef region is composed of a clay silt mixture. 
East of the inter-reef region the sediment is composed of 
fine to coarse grain sand particles [Ref. 13] . 

Based on the geo-acoustic modeling methods of Hamilton 
and Bachman [Ref. 14], a characterization of the reflective 
properties of the sediments can be accomplished. From Table 
1, the mean grain size can be determined from a description 
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of the sediment type. Grain size can then be correlated to 
porosity and saturated bulk density as illustrated in 
Figures 6 and 7. Porosity, which is a measure of the amount 
of voids occupied by water in a material, can then be used 
to determine the compressional wave speed of the sediment 
(Figure 8). The characteristic impedance, the product of 
the sediment density and wave speed, is an indicator of the 
reflective nature of a particular sediment composition. 

The geo-acoustic method described above indicates that 
the degree of reflectivity can be determined from the 
sediment type. In particular, reflectivity is proportional 
to grain size. It is anticipated, therefore, that the 
region on the eastern side of site 1 should be the most 
reflective and should have higher reverberation levels. The 
Howell Hook region should be less reflective, but still more 
reflective than the inter-reef region which should be the 
least reflective and have correspondingly lower 
reverberation levels . 

2 . Event Weather and Oceanography 

All acoustically relevant meteorological and 
oceanographic parameters were essentially invariant during 
the 12-hour duration of event 31, making for near- laboratory 
conditions. The sky was free of clouds with wind speeds 
less than 5 m/s. The sea surface was calm with a 
significant wave height of 0.5-0. 7 m. The surface current 
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measured during the event averaged around 0.25 m/s in a 
northerly direction. [Ref. 15] 

A total of 23 CTD observations were taken throughout 
site 1 between 6-9 February. Figure 9 shows the sound 
speed profiles (SSP) calculated from the CTD observations. 
Little spatial or temporal variability between the profiles 
is observed (< 1 m/s) , reflective of the near constant 
oceanic-atmospheric forcing. The mixed layer is seen to 
oscillate between 18 - 25 m. The downward refracting nature 
of the SSP ensures strong bottom interactions of the 
acoustic signal. [Ref. 15] 

3. Bottom Backscattering Measurements 

Bottom backscattering measurements taken with a 
vertical array during LWAD 99-1 [Ref. 16] indicated that the 
bottom scattering strength (S b ) exhibited a sin (9) 

dependence, not the sin‘(0) dependence of Lambert's Law. 
Previous results obtained from a shallow water area south of 
Long Island [Ref. 5] established a sin(0) bottom scattering 

strength relationship with a silty clay bottom having a 
surface sediment compressional wave speed of less than 1600 
m/s. Over sandy bottoms S b suggested a sin 2 (0) dependence in 

agreement with Lambert's Law. 

It is common practice to make backscattering 
measurements in just one location and then extend those 
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results to a wider area [Ref. 12], but this practice can be 
misleading. The spatial distribution of sediments in 
shallow water is generally highly variable over short 
spatial scale lengths. Scanlon et al . [Ref. 5] showed that 
the spatial scale can grade from hard sand to mud/silt clay 
sediments over distances less than 500 m. 

In this thesis, a bottom sediment map is constructed 
based upon spatial ly-varying changes in the beam 
reverberation data from an AN/SQS-53C tactical active sonar. 
While the results are considered preliminary due to a lack 
of high resolution bottom samples, they do indicate a high 
degree of spatial variability of sediment types that, only 
on the average, agree with the pre-exercise wide-area 
sediment assessment [Ref. 12]. It will be instructive to 
compare the results of the frequent sediment grab samples 
taken during LWAD 99-1 with the results of this thesis. A 
conclusion that can be stated now is: 

Single point bottom backscattering measurements can not 
predict the high spatial sediment variability and inverse 
techniques must be developed to estimate this variability. 

Another limitation of bottom backscattering analysis is 
the assumption that the scattering occurs at the sediment- 
water interface. Recent studies [Ref. 17] show that the 
scattering occurs from inhomogeneities in the sediment sub- 
bottom for slow (i.e., clay, mud) sediments at low grazing 
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angles (Figure 10). Therefore, a full wave model with an 
accurate three-dimensional model of the sub-bottom 
geoacoustic properties must be used to estimate TL in 
shallow water. Just a few years ago it did not seem 
possible to run a full wave model at 3 500 Hz, but with 
today's computing power, it is feasible. Within a few years 
it will be commonplace to run these full wave models in near 
real time and the real time IT sediment-mapping algorithm 
developed here is intended for this future application. 
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Sediment 

type 


No. 

samples 


Mean grain size 
mm phi size 


Sand, 

% 


Silt, 

% 


Clay, 

% 


Bulk 

grain 

density 

g/cm3 


Sand 
















Coarse 


2 


0.5285 


0.92 


100.0 


0.0 


0.0 


2.710 


Fine 


28 


0.1638 


2.61 


92.2 


4.1 


3.7 


2.709 


Very Fine 


16 


0.0988 


3.34 


81.0 


12.5 


6.5 


2.680 


Silty Sand 


40 


0.0529 


4.24 


57.0 


30.9 


12.1 


2.677 


Sandy Silt 


47 


0.0340 


4.88 


30.3 


57.8 


11.9 


2.664 


Silt 


19 


0.0237 


5.40 


7.8 


80.1 


12.1 


2.661 


Sand-silt-clay 


29 


0.0177 


5.82 


31.7 


42.9 


25.4 


2.689 


Clayey silt 


105 


0.0071 


7.13 


7.4 


58.3 


34.3 


2.656 


Silty clay 


54 


0.0022 


8.80 


3.9 


34.8 


61.3 


2.715 



Table 1. Continental terrace (shelf and slope) environment; 
Average sediment size analyses and bulk grain densities. 
From Ref . [14] . 
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II. TREATMENT OF DATA 



A. DESCRIPTION OF LWAD 99-1 REVERBERATION DATA 

The beam reverberation data used in this study was 
collected using the AN/SQS-53C fleet tactical sonar aboard 
the USS Vicksburg, during event 31 of exercise LWAD 99-1 on 
9 February 1999 [Ref. 10]. Beam reverberation data were 
analyzed from two modes of operation of the AN/SQS-53C, the 
SD (surface duct) and VD (variable depression) modes. For 
consistency, all pings that were analyzed used either a 
linear phase modulated pulse or a CP (coded pulse) waveform. 

The SD-mode uses a 300 millisecond pulse length and has 
a 360° azimuthal transmit coverage employing 72 receive 
beams. The VD-mode uses a 500 millisecond pulse length and 
has a 120° azimuthal transmit coverage employing 24 receive 

beams. In both modes, the 3-dB beamwidth is 5 degrees. The 
frequency range is from 2400 to 4500 Hz. The depression 
angle was 3°, while the range scale was set to 20 kyds . 

Event 31 was conducted such that the VD-mode pings were 
directed toward the center of the CG box. This was done to 
get complete coverage of the box and have RL information 
about the same geographic location from different 
directions . When the CG was traversing the western portion 
of the box, on a course of 000°T, the sonar center section 
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true bearing was 090°T (Figure 11) ; during the northern leg 
(course 090°T) the sonar center section true bearing was 

180°T (Figure 12), etc. The center section true bearing for 

the VD-mode is the bearing between beams 12 and 13 . 

All SD-mode pings (omni-directional) had a center 
section true bearing of 180°T. The center section true 
bearing for the SD-mode is the bearing between beams 36 and 
37. Beam 1 would therefore cover from 000°T to 005°T, beam 

2 from 005°T to 010°T, etc. Beam 72 would cover bearings 
355°T to 360°T. 

B. HYPOTHESIS 

Reverberation is a combination of backscattered energy 
from within the water column (volume) , the sea surface and 
the sea floor. For this study it is assumed that the volume 
and sea surface components of reverberation were constant. 
Therefore, it is hypothesized that for a single ping from 
the AN/SQS-53C different RL beam (bearing) time series 
should only differ because of variations in backscattered 
energy from differing bottom types. This hypothesis is 
further extended to encompass a set of pings (as many as 41 
pings were used in one data set) , as long as they occur in 
the same general area and there are no changes to the ping 
type. 
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The following assumptions were made to support this 
hypothesis : 

(1) The sound speed profile in the area around the 
source (within -20 kyds) is constant. 

• As seen in Figure 9, there is little spatial or 
temporal variability of the sound speed profile. 

(2) The output source level (SL) from the AN/SQS-53C in 
any direction does not vary (i.e., the AN/SQS-53C SL is 
independent in azimuth) . 

• The difference in performance of the AN/SQS-53C in 
any direction is considered negligible for this 
experiment. Analysis of individual RL time series, 
though, indicates that within 45° of the CG stern an 

increase in ambient (background) noise (AN) due to 
own ship noise is observed (Figures 13 and 14) . 
This phenomenon was accounted for in the processing 
of the RL data. 

(3) There are no major topographic changes that would 
significantly change the RL curve (i.e., no 
ups lope /downs lope effects, no pinnacles, etc.) 

• Based on the geologic characterization, there were 
no bottom topography features that would 
significantly effect the received reverberation 
levels . 
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C. PROCESSING METHOD 

An example of how the sediment effects reverberation 
can be seen by comparing two different RL time series in 
which distinctly different characteristics can be seen. 
Figure 13 is a single beam (beam 17) RL time series bearing 

082. 5°T. One sees that the RL curve reduces to the 

background noise level at around 18 s where the AN level at 
3.5 kHz is approximately 3 5 dB. In contrast, Figure 15 
(beam 37, bearing 182. 5°T) shows an example with the same AN 

level, for which the RL curve does not reach the AN level 
until about 23.5 s. Clearly the beam data is RL-limited 
longer in Figure 15 . It is suggested that this heightened 
RL is due to the presence of a more reflective bottom. The 
relationship between RL and AN as a function of detection 
range is illustrated in Urick [Ref. 3, Figure 2.2]. 

The general method proposed for RL curve analysis is as 
follows: (1) determine a characteristic RL curve for an area 
determined from the average of many pings and (2) compare 
individual beam data to determine its difference from the 
area-wide mean. 

A variety of statistical parameters derived from the 
beam RL time series can be examined to determine if the 
characteristics of the RL curve bear any correlation to 
sediment type. The statistical parameters investigated 
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include the mean slope, standard deviation, skew, and excess 
kurtosis . 



Before 


continuing. 


some background 


about 


the 


statistical 


method used in 


this thesis 


is required. 


The 


statistics 


calculated are 


not based on 


the 


entire 


time 


series but 


rather from 


segments of 


the 


time series 


partitioned 


into 0.5, 1 , 


or 2 second 


increments . 


In 



addition, a degree of overlap of the data is often useful, 
e.g., 0.25, 0.5, 0.75 (or 25/50/75% of the interval). The 
statistics derived from the overlapped segments are referred 
to as "moving" statistics (moving average, moving skew, 
etc.) . For example, a 25-second time series, with 320 data 
points per second, contains 8000 data points. With an 
interval/ overlap combination of 1/0.75 the first data point 
of the moving mean would be the mean of points 1 to 320, the 
second point would be the mean of points 81 to 400, the 
third would be the mean of points 161 to 480, etc. Figure 
16 shows an individual RL time series (black) overlaid by 
the 1/0.75 moving mean (white) RL curve. 

1. Example Using Slope Data, SD-Mode, 1/0.75 Moving 
Statistics 

From a data set of 41 pings an area-wide average slope 
is calculated by first taking the 1/0.75 moving mean for all 
of the beams (72) in each transmitted ping. The slope of 
each beam is determined by calculating the gradient (d/dt) 
of the mean RL curve. Calculating the mean from all (72x41) 
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slope data sets generates the area-wide average slope data. 
The area-wide average RL slope for all the pings examined is 
comprised of over 2500 individual RL slope time series. 
Figure 17 shows the area-wide average (solid line) RL slope 
and its standard deviation ( — ) . 

Figure 17 also shows the beam slope for a particular 
beam/ping relative to the area-wide average slope. Each 
point along the single beam slope has an associated range 
and bearing which corresponds to a geographic position. By- 
determining the difference between the area-wide average 
slope and the individual beam slope, a geographic contour 
plot can be generated to display the anomaly from the mean 
(Figure 18). Subsequently, different pings can be compared 
to determine if any geographic correlation exists between 
sets of pings to identify a unique but consistent area 
exhibiting similar RL characteristics. 

The other statistical parameters (STD, skew, and excess 
kurtosis) can be analyzed using an identical method to that 
used for the slope. 
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III. ANALYSIS 



A. ANALYSIS OF RL SLOPE DATA 

1 . SD-Mode 

A thorough analysis of all the RL slope plots was 
conducted to create a geographic representation of sediment 
type based on consistent spatial variations in the RL 
slopes. From a continuous sequence of 41 SD-mode pings a 
contour plot was created based on the deviations of the 
individual RL slopes from the area-wide average. 

Through the analysis of the slope data and comparisons 
to that from different pings, RL slope features can be seen 
that are consistent in their geographic location. The RL 
curves from ping 34 exhibit slope characteristics that are 
common in all pings in the set and will be discussed in 
detail. Figure 18 is a geographic contour plot of the 
deviations of each of the 72 azimuthal beams of ping 3 4 RL 
slope from the area-wide average RL slope. It is 
anticipated that regions where the deviations (from the 
average RL slope) are positive correspond to more reflective 
areas, while regions with a negative deviation are less 
reflective. For example, Figure 19 shows the RL and slope 
data of ping 34, beam 22. At approximately 8.5 s a region 
of increased RL is observed with a corresponding increase in 
RL slope. Figure 20 depicts the region of beam 22 enlarged 
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to show the RL slope details. The same geographic region of 
increased RL can be seen on contour plots of different pings 
(Figures 21 - 23). Note that the RL slope anomaly seen here 
is consistent in its geographic location indicating that it 
is due to a geographic feature, and not to a random 
occurrence . 

Figure 24 is similar to that of Figure 18 but is 
derived from a composite plot of all 41 SD-mode pings. It 
is the average deviation of RL slope from the area-wide 
average RL slope at each geographic location. Being an 
average of the RL slope deviations, this plot will show RL 
slope features that have geographic consistency, while 
eliminating those features that are not geographically 
consistent. Note that the area within + /- 45° of the CG 

stern were not included in the composite map and accounts 
for the notches, especially to the west, seen in the figure. 

An analysis of Figure 24 shows clear similarities to 
the pre-exercise sedimentary characterization (Figure 5). 
There are 3 distinct geographic locations observed in Figure 
24; east of 84.05°W is the f oraminif eral sands area, Howell 

Hook is west of 84.3°W and the middle is the inter-reef area 

which is interspersed with low and high reflective zones. 
The highly reflective zones could be uncovered (or thinly 
covered) reef segments whereas the less reflective spots may 
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be due to depressions/basins filled with unconsolidated 
sediments. It is notable that the foramini feral sands and 
Howell Hook regions are relatively uniform, while the inter- 
reef region shows high spatial variability of sediment 
types . 

2 . VD-Mode 

Due to the limited data available the VD-mode data were 
not used to create a composite map of the exercise area. 
However, individual pings were analyzed to see if there were 
any geographically consistent RL slope features as 
determined from the analysis of the SD-mode data. 

Due to the nature of the exercise the VD-mode pings 
were directed towards the center of the CG box. This was 
done to obtain full coverage of the area within the box but 
unfortunately created limitations for this study. In 
particular, sedimentary features such as Howell Hook and the 
foraminif eral sands region, which were beyond the boundaries 
of the exercise box, would not visible on the RL plots. 

The VD-mode pings did, however, demonstrate similar 
geographically consistent features as seen in the SD-mode. 
Figures 25 and 2 6 are individual VD-mode pings from which 
geographic features can be seen. Note that the high and low 
bottom reflection areas observed in Figures 25 and 26 are 
also seen at the same location in the SD-mode composite plot 
(Figure 24) . 
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B. HOS ANALYSIS 

The same analysis methods used on the RL slope data 
were employed to identify the correlation of HOS to sediment 
type. Figures 27-29 are plots of RL and the corresponding 
STD, skew, and excess kurtosis. The RL STD and skew plots 
show deviations that can be correlated to RL features. The 
RL kurtosis deviations, however, appeared random and showed 
no apparent correlation to sediment type. The observed RL 
STD and skew features are marginally useful in that they 
only provide an inference that some sedimentary change is 
present . They do not provide information about the 
characteristics of that feature (i.e., whether or not it is 
due to a more or less reflective sediment) . 

In order to further investigate RL statistics, a 
distribution analysis was conducted. Figure 30 is a plot of 
skew versus excess kurtosis for the detrended RL beam data 
of ping 189. The skew and excess kurtosis data shown are 
highly correlated with negative skew (characterized by an 
extended low amplitude tail compared to a Gaussian curve) 
correlated to positive excess kurtosis (characterized by a 
peaked distribution with more data centered around the mean 
compared to a Gaussian distribution) . Preliminary 
distribution analysis indicates that the reverberation data 
is not only non-Gaussian, but the data spans several 
distribution types [Ref. 18]. 
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It is important to note that statistical properties of 
the detrended beam time series are not commonly analyzed in 
classical detection theory applications. Instead, the 
output of a replica correlator is normally the input to a 
statistical analysis, and normalizations, such as whitened 
matched filters, etc., attempt to remove the amplitude 
decay, or slope, of the beam reverberation time series from 
the analysis [Ref. 18] . In this thesis we have taken a 
quite different approach by using the beam time series slope 
to invert bottom properties with the HOS of the detrended 
time series as derived product. If the thesis objective 
were focused on detection performance, we could have done a 
replica correlation analysis of the residual, or detrended 
beam time series. Instead, we are simply trying to 
correlate the HOS of the detrended time series to the 
environment. Since little is known about the HOS of the 
detrended beam time series, we are just reporting that they 
were highly non-Guassian and further investigation 
concerning their distribution function is needed. Although 
the HOS did not correlate well with sediment type, this does 
not imply that they will not correlate with some other 
environmental parameter, such as water depth or sound speed 
profile . 
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C. DISCUSSION 

1. Tactical Implications of Sediment Map 

The high variability of bottom sediment type in shallow 
water exerts a strong tactical significance because it can 
unexpectedly affect, often times adversely, the results of 
an acoustic search which is executed under the assumption 
that the bottom sediment is homogenous. Figure 31(A) is 
indicative of a search plan that was executed assuming that 
detection ranges were constant (regardless of sediment 
type) . This "mowing the lawn" technique could result in a 
lower probability of detection (PD) than anticipated due to 
regions where detection ranges are degraded because of a 
less reflective bottom. Prior to an acoustic search an 
environmental update could be conducted which includes an RL 
survey from which a sediment map, similar to the one in 
Figure 24, could be used to modify the search plan. 
"Searching by sediments", as seen in Figure 31(B), would 
account for detection range differences due to the bottom 
sediment type. In this case, the acoustic search would be 
executed under the appropriate, spatially- varying PD. 

2. Tactical Significance of Non-Gaussian HOS 

The analysis of the detrended beam reverberation data 
indicated that the backscattered energy is non-Gaussian, 
i.e., the values of skew and excess kurtosis are non-zero, 
and this also has important tactical implications. The 
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signal processing design incorporated into the AN/SQS-53C 
and other sonar systems is based on the assumption that the 
RL distribution is Gaussian. Challenging the Gaussian 
assumption raises two important questions: 

(1) Is the non-Gaussian nature of the beam 
reverberation HOS due to environmental effects (e.g., 
shallow water acoustic propagation) or because of sonar 
signal processing design? 

(2) If the RL HOS are due to the environment, will the 
HOS change when calculated for a different environment 
(e.g., deep water, continental slope or soft, muddy 
bottoms) ? 

The first question could be answered by performing the 
same statistical analysis of AN/SQS-53C RL data but taken 
from another environment. If the HOS are similar in varied 
environments, the nature of the HOS is probably due to the 
signal processing equipment. In this case, the non-Gaussian 
issue would be resolved by making equipment modifications. 

If the non- zero HOS are in fact caused by environmental 
factors, then the existing signal processing algorithms, 
from which sonar system design is based, are not valid. 
Currently correlation detectors or . energy threshold 
detectors are optimized, in the probability of 
detection/probability of false alarm (PD/PFA) sense for 
Gaussian beam reverberation statistics. If the signal is 
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not Gaussian, then non-Gaussian treatment of the RL data 
must be used to process the reverberation. Solinsky [Refs. 
19-21] has shown that non-zero HOS algorithms can be 
designed which will efficiently mitigate the effects of 
reverberation . 

A fundamental conclusion from this study is that the 
statistical distribution and its associated properties of 
the reverberation HOS needs to be determined in-situ and 
incorporated into the PD/PFA algorithms to improve echo to 
reverberation processing of the AN/SQS-53C sonar system. 
The HOS characteristics of the environment can be quantified 
through analysis of the abundant AN/SQS-53C data available 
in the Naval Undersea Warfare Center (NUWC) active sonar 
data base. 

The non-Gaussian reverberation statistics problem can 
be addressed in several ways, including the development of a 
neural net (NN) algorithm, schematically illustrated in 
Figure 32. The NN approach can be adapted to the changing 
environment and this technique is highly recommended for 
further development for the AN/SQS-53C. 

Another deterministic approach to addressing non- 
Gaussian statistics is to apply matched field processing 
(MFP) concepts. Here the transmitted signal replica is 
propagated through the shallow water environment using an 
accurate propagation model based upon accurate environmental 
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parameters determined by ITs such as the one developed in 
this thesis. Any procedure to process the transmitted 
waveform to accurately reflect the non-Gaussian 
characteristics of the environment will improve the 
correlation of the received energy which has been propagated 
in the environment and reflected off the target and other 
scatterers . If the transmitted waveform is not similar to 
that of the received energy, the correlation will be low and 
detection performance poor. On the other hand, if the 
volume, bottom and surface scattered returns correlate "too 
well" with the transmitted waveform, the active sonar system 
will be overloaded with false targets. MFP or neural net 
algorithms are ways to address the problem of non-Gaussian 
statistics by including the realities of propagation and 
scattering in shallow water in the processing of the 
transmitted waveform prior to correlation with received beam 
data. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



Exercise LWAD 99-1, which took place in February 1999, 
included an inversion technique experiment during which 
AN/SQS-53C reverberation data was collected. The 
reverberation data was used to investigate the possibility 
of employing inversion techniques to infer bottom 
sedimentary properties. 

The hypothesis presented in this thesis was based on 
the assumption that the magnitude of backscattered energy 
and its rate of spatial decay were directly related to the 
acoustic reflectivity of the seabed. By examining the slope 
of the RL curves from various sedimentary provinces, one 
should be able to identify and associate various reverberant 
areas with unique sediment types. This study examined 
various reverberation level statistics including RL slope, 
standard deviation, skew, and excess kurtosis, to test this 
hypothesis. Of the various statistical parameters the RL 
slope was found to be most clearly aligned with the degree 
of bottom reflectivity. Over less reflective (mud/silt) 
bottoms a decrease in RL slope was noted while acoustic 
interactions with a more reflective bottom (sand) causes an 
increase in RL slope. Based on this finding, a technique 
was developed using the deviation of the RL slope for a 
particular ping/beam from an area-wide average RL slope to 
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generate a geographic map illustrating bottom sedimentary 
characteristics . 

Analysis of higher order statistics (HOS) , which 
included standard deviation, skew and excess kurtosis, 
revealed no clear correlation to sediment type. Variations 
in standard deviation and skew did provide indications of 
sedimentary change . 

An analysis of the statistical distribution of the HOS 
indicated that the RL data did not derive from a Gaussian 
population. It also showed that the distribution did not 
appear to come from any well-known distribution (e.g., 
Poisson, log normal, Rayleigh, etc.). The non-Gaussian 
nature of the RL data has serious implications for sonar 
signal processing because most signal processing algorithms 
assume a Gaussian distribution to the noise field. 

It is recommended that the RL slope technique for 
generating a sediment map be used by fleet assets as part of 
efforts to develop in-situ and archived bottom reflection 
databases and that these results be incorporated into 
shallow water ASW search planning. In addition, further 
investigation into the nature of RL statistics and its 
impact on sonar system performance should be undertaken. 
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APPENDIX A. FIGURES 




TimeRange 



Surface scattering 




Figure 1. Schematic and characteristic RL curve for a 
deep water environment. Reverberation is due primarily 
to backscattering from volume and sea surface scatterers . 
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Time/Range 



Surface scattering 




Shallow Water 



Figure 2 . Schematic and characteristic RL curve for a 
shallow water environment. Reverberation is generally 
louder being a function of volume, sea surface, and 
bottom backscatterers . High RL levels often mask the 
target echo. 
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Figure 3. LWAD 99-1 Exercise Area. Event 31 (Inversion 
Techniques) was conducted at site 1. Depths in meters. 
From Ref . [ 9 ] . 
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Figure 4. Event 31 experiment geometry. Over 900 pings 
of the AN/SQS-53C tactical sonar were completed between 
073 0Z and 1800Z on 9 February 1999. From Ref. [10]. 
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Figure 5. Exercise LWAD 99-1 surface sediment 
distribution. After Ref. [11]. 
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Figure 6. Mean grain size versus porosity. Round dots 
represent continental terrace (shelf and slope) samples. 
From Ref . [ 14 ] . 
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MEAN GRAIN SIZE 



Figure 7. Mean grain size versus saturated bulk density. 
Round dots represent continental terrace (shelf and 
slope) samples. From Ref. [14]. 



39 



1900 



o 

o 

_1 

LU 

> 

Q 

Z 

D 

O 

co 



1860 
1820 
1780 
1740 
1700 
1660 
1 620 U 

1580|> 

1540 

1500h 



• • 

V- 

+ * 

< • 
- «• • 



I*. 



• *<; - 

. • S' • • 

> v • • 9 
• 

• •• • 

• * • • 

•• • 

•-V . V 






I I I I 



■ ’■ u 



30 40 50 60 70 80 90 100 

POROSITY. % 



Figure 8. Porosity versus sound velocity, continental 
terrace (shelf and slope) . From Ref. [14] . 
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SITE 1 CTD DOWNCASTS 



Sound Speed (m/sec) 




Figure 9. 23 total CTD downcasts were taken at site 1 
from 2/6/99 to 2/9/99. Bottom interaction of acoustic 
energy is ensured due to downward refracting nature of 
sound speed profile (SSP) . From Ref. [15] . 
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Figure 10. Scattering mechanisms are from (1) the water- 
sediment interface, (2) sub-bottom volume, and (3) the 
sediment-basement interface. From Ref. [17], 
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Figure 11. VD-mode ping with center section true bearing 
of 090°T. VD-mode pings are steered toward the center of 
the CG box. 
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Figure 12 . VD-mode ping with center section true bearing 
of 180°T. VD-mode pings are steered toward the center of 
the CG box. 
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Ping: 189, Mode: SD. Beam: 17, Brng(T): 82.5 
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Figure 13. RL time series for a forward projected beam. 
RL reduces to the ambient noise (AN) level of 
approximately 35 dB at about 18s. Compare with Figure 14 
(stern beam) to see difference in AN level. 
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Ping: 189, Mode: SD, Beam: 54, Brng(T): 267.5 



® 100 



<r» 

> 

<D 



C 

o 



<D 

-2 

O) 

> 

CD 

or 




10 15 

Time (sec) 



Figure 14. RL time series of stern beam. The ambient 
noise (AN) level is approximately 20 dB greater than non- 
stern beam due to ship noise. 
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Ping. 189, Mode: SD, Beam: 37, Brng(T): 182.5 




Figure 15 . RL curve reduces to background noise level at 
around 23.5 s. Compared to figure 13 the longer RL decay 
time is associated with backscattering from a more 
reflective bottom. 
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Reverberation Level (dB) 



Ping: 189, Mode: SD, Beam: 17, BrngfT): 82.5, Interval: 1, Overlap: 0.75 




Figure 16. Individual RL time series (black) overlaid by 
the 1/0.75 moving mean (white) RL curve. 
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Ping: 189, Mode: SD, Beam: 17, Brng(7): 82.5, Interval: 1 , Overlap: 0.75 
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Figure 17. Slope of the 1/0.75 moving mean (-*) for an 
individual beam compared with the area-wide average slope 
(solid line) and its standard deviation (--) . 
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Ping:0034 .Mode: SD. Interval: 1, Overlap: 0.75 
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Figure 18. Ping 34. Geographic contour plot of 
deviations of all 72 beams of slope data from the area- 
wide average slope. Warmer colors represent more 
reflective sedimentary areas than the cooler colors. 
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Figure 19. Ping 34, beam 22. An increase in RL at 8 . 5 s 
corresponds with an increase in RL slope. 
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Figure 20. Ping 34, beam 22. Enlarged contour plot to 
display the region of increased RL. 
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Ping: 25 , Mode. SD, Interval: 1 , Overlap: 0.75 
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Figure 21. Ping 25. Note region of increased RL at same 
geographic location as Ping 34. 
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Figure 22. Ping 30. Note region of increased RL at same 
geographic location as Ping 34. 
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Figure 23. Ping 35. Note region of increased RL at same 
geographic location as Ping 34. 
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Deviation from average RL slope 




Figure 24. Composite plot of deviation of RL slope from 
the area-wide average for 41 SD mode pings. Warmer 
colors represent more reflective sedimentary areas than 
the cooler colors. 
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Figure 25. VD-mode ping 2. Note that the features 
pointed out are also seen in the SD-mode composite plot. 
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Figure 26. VD-mode ping 4. Note that the features 
pointed out are also seen in the SD-mode composite plot. 
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Figure 27. Ping 34, beam 22. An increase in RL at 8 . 5 s 
corresponds to an increase in RL STD. 
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Figure 28. Ping 34, beam 22. An increase in RL at 8.5 s 
corresponds to an RL skew feature . 
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Ping: 034, Mode: SD, Beam: 23, Brng(T): 292.5, Interval: 1 , Overlap: 0.75 




Figure 29. Ping 34, beam 22. The RL kurtosis features 
appear random with no apparent correlation to RL. 
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Figure 30. Skew vs. excess kurtosis for the detrended 
beam reverberation data of ping 189. 
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Figure 31. Old search method (A) vs. new search method 
(B) with knowledge of local sediment types. 
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From Ref. 



[ 21 ] . 
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APPENDIX B. MATLAB FILES 



The matlab program files enclosed are for analyzing 
beamform data from exercise LWAD 99-1. The program files do 
not include the actual beamform ping data, or the ping 
parameter files. The matlab files included here are for 
guidance in creating a graphical user interface (GUI) for 
analyzing beamform data. 

Beamform data files contain the following variables: 
bandwidth - bandwidth (Hz ) 
bmf - complex beamform data 
centerfreq - center frequency (Hz) 
ping_date - date (e.g. '02/09.99') 

ping_time - time (e.g. '11:18:15') 

pulselength - pulse length (ms) 
samplerate - samplerate (Hz) 

The parameter file 'evt31.mat' contains a vector of 
times and corresponding CG and TGT positions from LWAD 99-1 
event 3 1 : 

evt31_time - time in seconds from midnight 
evt31_date - event 31 date (e.g. '02/09.99') 

cg_lat - vector of CG latitudes 
cg_lon - vector of CG longitudes 
tgt_lat - vector of TGT latitudes 
tgt_lon - vector of TGT longitudes 
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The parameter file ' evt31prm . mat ' contains a vector of 



times and corresponding ping parameters : 
botdepth - bottom depth (fm) 
course - CG course (°T) 

ctrbrng - Center of ping section bearing (°T) 

deangle - sonar depression angle 
speed - CG speed (kts) 

timeday - time in seconds from midnight 
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%LWAD 

% 



% 

% 

% 

% 

% 

% 

% 

% 



% 

% 
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% 

% 

% 

% 



% 

% 



Run a GUI to analyze LWAD 99-1 data 

A. Begin by selecting a ping to analyze 
The following will occur: 

(1) A map of the LWAD 99-1 exercise will be created 
showing the following: 

(a) blue eg track history, red circle is eg 
current location 

(b) red tgt track history, blue circle is tgt 
current location 

(c) ping section is shown in yellow with dashed 
yellow center 

section true bearing line. Yellow asterisks 
show lOnm range. 

(d) latitude/ longitude grid lines (black dashed 
lines ) 

(e) Sediment map showing Howell hook (black circles) 
with phi size approx 4-7, clay silt (black dots) 
with phi size approx 6-8, and sand 

(black slashed lines) with phi size approx 1. 
This map is based on data received from 
J. Fulford (NRL Stennis) . 



% (2) A list of ping parameters will be displayed in 

% the 'LWAD Parameters' section. 

% 

% (3) A data analysis section will be created 

% in the GUI. 



% B. To create a contour plot press 'Contour Plot' button. 
% - This plot is created by determining the 

% difference between area averaged slope data and 

% slope of individual beam data. 

% - The plot is a filled contour plot. 

% - A regular (non-filled) contour plot may be 

% created by using the function ' lwad_contour .m' . 

% 

% Note: when the 'contour plot' button is pressed, the 
% contour plot is created 

% using an interval of 1 second and an overlap of 

% .75. The interval /overlap is not what is selected 

% in the LWAD Data Processing Program GUI. If you 

% want to investigate different interval /overlap 

% combinations, you'll probably want to run 

% ' lwad_contourf ' function separate from the LWAD 
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Data Processing Program GUI. 



% 

% 

% 

% 

% 



% 

% 

% 

% 



% 

% 

% 

% 

% 

% 

% 

% 

% 



C. To analyze beamform data 

(1) select the type of plot(s) to display 

- 201og|data| - is RL data in dB with averaged 
data in red 

- Slope - slope of RL data in blue with area 

averaged slope in red. Dashed red 
line is +/- 1 standard dev. 

- STD - standard deviation in blue with area 

averaged std in red. Dashed red line 
is +/- 1 standard dev. 

- Skew - Skewness in blue with area averaged 

skewness in red. Dashed red line is 
+ /- 1 standard dev. 

- Kurtosis - Excess Kurtosis in blue with area 

averaged excess kurtosis in red. 
Dashed line is + /- 1 std. 

(2) Select the beam you would like to analyze 
1-24 (VD) or 1-72 (SD) 

(3) Select the interval and overlap to be used by 

mmean.m, mstd.m, mskewness.m, mkurtosis.m. I 
usually use 1 sec and .75 overlap. Other 
combinations may be used to see the effect. 

(4) Press 'Finish' to see the plot(s). 

(5) Press 'Close plots' to close all plots but the 

event map, and the LWAD GUI. 



% D. To analyze another file press 'Browse'. 

% 

% E. Press 'Quit' to end your session and close all 
% open plots . 

% 

% See also LWAD_CONTOUR , LWAD_CONTOURF , LWAD_MAP, 

% LWAD_POS IT , LWAD_PRM , MMEAN, MSTD, MKURTOSIS, 

% MSKEWNESS . 



% Created: May 8, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 



menu_num = 1 ; 
figure 

load 'e:\schalm\lwad 99-l\data_drive ' 
lwad_menul 
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% LWAD_MENU 1 



% 

% This will create a UI to select a file for the LWAD 
% program 
% 

% See also LWAD 

% Created: April 29, 1999 

% David Schalm 

% 

% Last Modified: July 26 , 1999 
% David Schalm 



fig = gcf ; 



set (fig, . . . 

'Units' , 'points' , ... 

'Color', [0.8 0.8 0.8], ... 

'MenuBar' , 'none' , ... 

' Number Title ' , 'off', ... 

'Name' , 'LWAD Data Processing Program' , ... 

' Position' , [450 240 300 300], ... 

' Resize ' , ' off ' ) ; 
uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ' , ' lef t ' , ... 

' Position' , [10 280 180 10], ... 

' String' ,' LWAD Data file', ... 

'FontWeight' , 'bold' , ... 

' Style ' , ' text ' ) ; 

hO = uicontrol ( 'Units ', 'points ' , ... 

'BackgroundColor' ,[111], ... 

' HorizontalAlignment ' , ' left ' , ... 

'Position' , [10 260 180 15] , ... 

' String ', [da ta_drive '*.mat'], ... 

' Style' , ' edit ' ) ; 

callbackl = [' [data_f ile , data_drive] = ' . . . 

' uigetf ile (data_drive, ' ');' ... 

' set ( fig, ' ' ' ' Name ' ' ' ' , ' ' ' ' Loading File . . . ' ' ' 

');' 'if data_file == 0, load ' ... 

' ' 'e:\schalm\lwad 99-l\data_drive ' ' ' ';'... 
'set(fig, ' '''Name''' ',' ... 

' ' 'LWAD Data Processing Program' ' ' ' ) ; ' ... 

' , else, lwad_menu2 , end' ] ; 
uicontrol ( 'Units ', 'points ' , ... 

'Callback' , callbackl, . . . 

' Position' , [200 260 40 15], ... 
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'String' , 'Browse' ) ; 

callback2 = [ ' ButtonName=questdlg ( ' ... 

' ' 'Are you sure you want to quit?' ' ' . . . 

',' ' ''Quit LWAD' ' ' ',' '''Yes''' ',' ... 

' ' 'No' ' ' ' , ' ' ' 'Yes' ' ' ');'... 

'if length (Butt onName) == 3, close all, end']; 
uicontrol ( 'Units ', 'points ' , ... 

' Callback' , callback2 , ... 

' Position' , [10 10 40 15], ... 

' String' , 'Quit' ) ; 
uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ' , [0 . 5 0.5 0.5], ... 

' HorizontalAlignment ' , ' lef t ' , ... 

' Position' , [5 250 290 1], ... 

'String',' ', ... 

' Style ' , ' text ' ) ; 
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%LWAD_MENU2 



% 

% This will display the file parameters, and add a UI 

% section for selecting plot(s) and plot options. 

% 

% Note: when the 'contour plot' button is pressed, the 
% contour plot is created using an interval of 1 

% second and an overlap of .75. The interval/overlap 

% is not what is selected in the LWAD Data Processing 

% Program GUI. If you want to investigate different 

% interval /overlap combinations, see 

% ' lwad_contourf .m' . You'll probably want to run 

% ' lwad_contourf ' separate from the LWAD GUI. 

% 

% See also LWAD, LWAD_C ONTOURF 

% Created: April 29, 1999 

% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 



fig = gcf ; 



save 'e:\schalm\lwad 99-l\data_drive .mat ' data_drive 
fname = [data_drive data_f ile] ; 
load ( fname , ' -mat ' ) 

[course, speed, ctrTbrng,botdepth, deangle] = lwad_prm ( fname ) ; 
[cglat , cglon, tgtlat , tgtlon] = lwad_posit ( fname) ; 
num_of_beams = size(bmf,2); 
set (hO , ' string ' , data_f ile) ; 

uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ' , ' lef t ' , ... 

' Position ',[ 10 230 180 10], ... 

'String' , 'LWAD Parameters' , ... 

'FontWeight' , 'bold' , ... 

' Style ' , ' text ' ) ; 
uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position ',[ 10 210 140 10], ... 

' String ',[' Ping date: ' ping_date] , ... 

' Style' , ' text' ) ; 
uicontrol (' Units ',' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 
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'HorizontalAlignment ' , ' left' , ... 

'Position' , [150 210 140 10], ... 

' String' ,[' Ping Time: ' ping_time] , ... 

' Style' , ' text' ) ; 

lwad_map ( f name ) ; 
figure ( fig) 

uicontrol (' Units ', 'points ' , ... 

' BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment ' , ' lef t ' , ... 

' Position' , [10 200 140 10], ... 

' String' ,[' CG Position: ' num2str (cglat ) 'N ' .. 

num2str (cglon) 'W'], ... 

' Style ' , ' text ' ) ; 
uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment' , 'left' , ... 

' Position' ,[ 150 200 140 10], ... 

' String' ,[' TGT Position: ' num2str (tgtlat) 'N 

num2str ( tgtlon) 'W'], ... 

' Style ' , ' text ' ) ,- 
uicontrol (' Uni ts ', 'points ' , ... 

' BackgroundColor ', [0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position' ,[ 10 190 140 10], ... 

' String' ,[ 'CG Course (T) /Speed (kts ) : ' ... 

num2str (course) '/' num2str (speed) ] , ... 

' Style ' , ' text ' ) ; 
uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment' , 'left' , ... 

' Position' , [150 190 140 10], ... 

' String' ,[' Center Bearing: ' num2str (ctrTbrng) ] , 

' Style' , ' text' ) ; 
uicontrol ( ' Units ' , ' points ' , ... 

'BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment' , ' left' , ... 

' Position' , [10 180 140 10], ... 

' String' ,[' Bottom Depth(fm): ' num2str (bo t depth ) ] 

' Style' , ' text ' ) ; 
uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment ',' left ' , ... 

' Position' , [150 180 140 10], ... 

' String' ,[' Depression Angle: ' num2str (deangle) ] , 

' Style' , ' text' ) ; 
uicontrol ( ' Units ' , ' points ' , ... 



72 



' BackgroundColor ',[0.8 0.8 0.8], ... 

' HorizontalAlignment ' , ' left' , ... 

' Position' , [10 170 140 10], ... 

' String' ,[ 'Bandwidth: ' num2str (bandwidth) ' Hz'], . 

' Style ' , ' text ' ) ; 

uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ' , ' lef t ' , ... 

' Position' , [150 170 140 10], ... 

' String' ,[' Center Freq: ' num2str (centerfreq) ... 

' Hz'],' Style ' , ' text ' ) ; 
uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position ',[ 10 160 140 10], ... 

' String' ,[' Pulse length: ' num2str (pulselength) ... 

' msec ' ] , ' Style ' , ' text ' ) ; 
uicontrol ( 'Units ', 'points ' , ... 

'BackgroundColor ', [0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position ',[ 150 160 140 10], ... 

' String ',[' Sample rate: ' num2str ( samplerate) ... 

' Hz' ] , 'Style' , 'text' ) ; 

% Create time vector 
data_start = 1; 

data_duration = length (bmf ) /samplerate - ( 1 /samplerate) 
+ data_start; 

time = (data_start:l/samplerate:data_duration)'; 

uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ',[0.8 0.8 0.8], ... 

'HorizontalAlignment' , 'left' , ... 

' Position ',[ 10 150 140 10], ... 

' String ',[' Data Duration: ' num2str (data_duration) . . 

' Sec ' ] , ... 

' Style ' , ' text ' ) ; 

uicontrol ( ' Units ' , ' points ' , ... 

'BackgroundColor' , [0 . 5 0.5 0.5], ... 

'HorizontalAlignment' , 'left' , ... 

'Position' , [5 140 290 1] , ... 

' String' , ' ' , ... 

' Style ' , ' text ' ) ; 

uicontrol (' Units ',' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 
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'Position' , [10 120 180 10], ... 

' String' , . . . 

'Select the plots you would like to display.', 
'FontWeight' , 'bold' , ... 

' Style' , ' text' ) ; 

hi = uicontrol ( 'Units' , 'points' , ... 

' BackgroundColor ' , [ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ' , ' lef t ' , ... 

' Position' , [10 100 60 10], ... 

'String', '20 log |data|', ... 

'Style', 'checkbox'); 
h2 = uicontrol ( 'Units' , 'points' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position' ,[ 10 85 60 10], ... 

' String ' , ' Slope ' , ... 

'Style', 'checkbox'); 
h3 = uicontrol (' Units ',' points ' , ... 

' BackgroundColor ', [ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment' , ' left' , ... 

' Position' , [10 70 60 10], ... 

'String' , 'STD' , ... 

'Style', 'checkbox'); 
h4 = uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position' , [10 40 60 10], ... 

' String ' , ' Kur tosi s ' , ... 

'Style', 'checkbox'); 
h5 = uicontrol ( 'Units ', 'points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

' HorizontalAlignment ',' lef t ' , ... 

' Position' , [10 55 60 10], ... 

' String' ,' Skew' , ... 

'Style', 'checkbox'); 
uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[ 0 . 8 0.8 0.8], ... 

'HorizontalAlignment ',' left ' , ... 

' Position' , [85 100 65 10], ... 

' String ',' Beam Selection:', ... 

' Style ' , ' text ' ) ; 
beam = {1:1: num_of_beams } ; 
h6 = uicontrol (' Units ',' points ' , ... 

'HorizontalAlignment' ,' left ' , ... 

'BackgroundColor' ,[111], ... 

' Position' , [85 40 65 55], ... 

' String ', beam, ... 
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' Style ' , ' listbox ' , ... 

'ListBoxTop' , 1); 

uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ', [0 . 8 0.8 0.8], ... 

'HorizontalAlignment ' , ' left ' , ... 

' Position' , [160 100 65 10], ... 

' String' ,' Interval (sec)', ... 

' Style' , ' text' ) ; 

interval_sel = {.5 1 2); 
h7 = uicontrol ( 'Units ', 'points ' , ... 

'BackgroundColor ', [1 11], ... 

'HorizontalAlignment' , 'left' , ... 

' Position' , [160 40 60 55], ... 

'String' , interval_sel, ... 

'Style', 'listbox'); 
uicontrol ( ' Units ' , ' points ' , ... 

' BackgroundColor ',[0.8 0.8 0.8], ... 

' HorizontalAlignment ',' left ' , ... 

' Position' , [230 100 65 10], ... 

'String' , 'Overlap (0-1)', ... 

' Style ' , ' text ' ) ; 

overlap_sel = {0 .25 .5 .75}; 
h8 = uicontrol ( 'Units ', 'points ' , ... 

'BackgroundColor' ,[111], ... 

' HorizontalAlignment ',' left ' , ... 

' Position' , [230 40 60 55], ... 

' String' , over lap_sel , ... 

'Style', 'listbox'); 

callback4 = ['set (fig,' '''Name''' ',' ... 

' ' 'Getting Contour Plot. . . ' ' ' ');'... 

' lwad_contour f ( f name , 1 , . 7 5 ) ; ' . . . 
'set(fig, ' '''Name''' ',' ... 

' ' 'LWAD Data Processing Program' ' ' ' ) ; ' ] 

uicontrol ( ' Units ' , ' points ' , ... 

'Callback', callback4, ... 

' Position' , [80 10 60 15], ... 

' String' ,' Contour Plot'); 
callback3 = ['set (fig,' '''Name''' ',' ... 

' ' ' Getting Plots ...''' ');'... 

' lwad_menu3 ' ] ; 

uicontrol ( ' Units ' , ' points ' , ... 

'Callback', callback3 , ... 

' Position' , [250 10 40 15], ... 

' String' , ' Finish' ) ; 

set(fig, 'Name' , 'LWAD Data Processing Program' ) 
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% LWAD_MENU 3 



% 

% This script file is used to create the selected plots 

% 

% See also LWAD 

% Created: April 29, 1999 

% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 

beam_s elect ion = beam{l} (get (h6 , 'value ')) ; 
data_file = get (hO ,' string' ) ; 
interval = interval_sel {get (h7 , 'value') }; 
overlap = over lap_sel {get (h8, 'value' )} ; 
spacing = interval - interval *over lap; 

N = round ( interval *samplerate) ; 

% Create time vector 
[mtime,ao] = mmean ( time, N, overlap) ; 
if ao ~= overlap 
overlap = ao; 

disp ( [ 'Note : actual overlap is ' num2str (ao) ] ) ; 

end 



%%%%%%% PLOT THE DATA %%%%%%%%%%%%% 

% Display plots of data/mean, slope, kurtosis, skewness 

if data_f ile ( 5 : 6 ) == 'SD' 

area_stats_fname = ['e:\schalm\lwad 99-l\area_stats\ ' . 
data_f ile ( 5 : 6 ) '_stats' ... 
num2 s t r ( interval ) ' .mat'] ; 

else 

area_s tat s_f name = ['e:\schalm\lwad 99-l\area_stats\ ' . 
data_file (19 : 20) '_stats' ... 
num2 s t r ( interval ) ' . mat ' ] ; 

end 

load ( area_s tat s_f name ) 

if get (h8 , 'value ' ) == 1 
area_stats = stats_0; 
elseif get (h8 , 'value ' ) == 2 
area_stats = stats_25; 
elseif get (h8 , 'value ' ) == 3 
area_stats = stats_5; 
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4 



elseif get (h8 value ' ) == 
area_stats = stats_75; 

end 

log_data = 20*logl0 (abs (bmf ( : , beam_s elect ion) ) ) ; 
thetal = ctrTbrng - num_o f _beams * 5 / 2 + 5/2; 
theta = thetal + 5* (beam_s elect ion - 1) ; 
if theta < 0, theta = theta + 360; end 

title_string = ['Ping: ' data_f ile (14 : 16) Mode: SD' ... 

' , Beam: ' num2str (beam_selection) . . . 

', Brng(T): ' num2str (theta) ... 

Interval: ' num2str (interval) Overlap: ' ... 

num2str (overlap) ] ; 

% plot 201og|data| 
if get (hi value ' ) == 1 
figure, b = gcf; 

plot ( time, log_data, ' k' ) , title ( title_string) , 
ylabel ( 'Reverberation Level (dB) ' ) , xlabel ( 'Time (sec) ' ) 
axis([0 data_duration 0 150]), grid on, hold on 
[mean_data] = mmean(log_data,N, overlap) ; 
plot (mtime, mean_data, ' w' ) 

end 

% plot the slope of 201og|data| 
if get (h2 , 'value ' ) == 1 
figure, b = gcf; 

[mean_data] = mmean(log_data,N, overlap) ; 
slope_data = gradient (mean_data) ; 

plot (mtime, slope_data/ spacing, ' k*- ' ) , title ( title_string) 
ylabel ('RL Slope '), xlabel (' Time (sec)'), grid on, hold on 
plot (mtime, area_stats (1, : ) /spacing, 'k' ) 
plot (mtime, (area_stats ( 1 , : ) -area_stats (2 ,:))/.. . 
spacing, ' --k' ) 

plot (mtime, (area_stats (1, : ) +area_stats (2, :))/.. . 
spacing, ' --k' ) 

legend ( [ 'Beam ' num2str (beam_selection) . . . 

' Slope' ] , 'Average Slope' , 'Average Slope STD' ) 

end 

% detrend the data w/ a 3rd order polynomial 

log_data = log_data-polyval (polyfit (time, log_data, 3) , time) ; 

% Plot the STD 
if get (h3 ,' value ' ) == 1 
figure, b = gcf; 



77 



[std_data] = ms td ( log_data, N, overlap) ; 
plot (mtime, std_data, 'k*-' ) , title ( title_string) 
ylabel ( ' STD' ) / xlabel ( ' Time (sec)') 
grid on, hold on 
plot (mtime, area_stats (3 , : ) , 'k' ) 

plot (mtime, (area_stats (3 , : ) -area_stats (4 , : ) ) , ' --k' ) 
plot (mtime , (area_stats (3 , : ) +area_stats (4 ,:)),' --k' ) 
legend ([' Beam ' num2str (beam_selection) ' STD'],... 
'Average STD', 'STD of Average STD') 

end 

% Plot the Kurtosis 
if get (h4 , 'value ' ) == 1 
figure, b = gcf; 

[kurtosis_data] = mkurtosis (log_data, 1,N, overlap) ; 
plot (mtime, kurtosis_data, ' k*- ' ) , title ( title_string) 
ylabel (' Kurtosis ') , xlabel ('Time (sec)') 
grid on, hold on 

plot (mtime, area_stats (7 , : ) , ' k' ) 

plot (mtime, (area_stats (7 , : ) -area_stats (8 , : ) ) , ' --k' ) 
plot (mtime , (area_stats ( 7 , : ) +area_stats (8 , : ) ) , ' --k' ) 
legend ([' Beam ' num2str (beam_selection) ' Kurtosis'],. 
'Average Kurtosis ',' STD of Average Kurtosis') 

end 

% Plot the Skew 
if get (h5 , 'value ' ) == 1 
figure, b = gcf; 

[skew_data] = mskewness (log_data,N, overlap) ; 
plot (mtime, skew_data, 'k*- ' ) , title ( title_string) 
ylabel ( ' Skew' ) , xlabel ( ' Time (sec) ' ) 
grid on, hold on 

plot (mtime, area_stats ( 5 , : ) , ' k' ) 

plot (mtime , (area_stats ( 5 , : ) -area_stats (6,:)),' — k') 
plot (mtime , (area_stats (5 , : ) +area_stats ( 6 ,:)),' --k' ) 
legend ([' Beam ' num2str (beam_selection) ' Skew'],... 
'Average Skew', 'STD of Average Skew') 

end 

set ( fig, ' Name ',' LWAD Data Processing Program') 
figure (fig) 

uicontrol ( 'Units ', 'points ' , ... 

'Callback', [ ' closebut ( [ f ig 100])'], ... 

' Position' , [165 10 60 15], ... 

' String ',' Close Plots'); 
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function [area_stats] = lwad_s tats (file_type, interval .. . 

, overlap) 

' %LWAD_STATS Determine the statistics for the LWAD 99-1 
% exercise area. 

% [AREA_STATS,ERROR_LOG] = LWAD_S TATS ( F I LE_TY PE , 

% INTERVAL, OVERLAP) will return the statistics for the 
% selected parameters. The function will calculate 
% statistics for each individual beam, average these 
% statistics, and determine the standard deviation. 

% 

% FILE_TYPE is either 'sd', 'cp', or ' cw ' ... depending 

% on which mode you want to analyze . 

% 

% INTERVAL and OVERLAP are required by funtions MMEAN, 

% MSKEWNESS, and MKURTOSIS to calculate the respective 
% statistics. If no INTERVAL or OVERLAP are specified, 

% the program defaults to INTERVAL = 1, OVERLAP = .75. 

% 

% The mean and std for the following statistics are 
% determined: 

% (1) Slope 

% ( 2 ) STD 

% (3) Skewness 

% (4) Excess Kurtosis 

% 

% Note: STD, Skewness, and Excess Kurtosis are determined 
% from 'detrended' data. The data is detrended by 

% subtracting a best fit 3rd order polynomial from 

% the data. Slope is not detrended. 

% 



% 


The 


output 


file AREA_STATS is organized as such: 


% 




Row 


1: 


mean of the slope 


% 




Row 


2 : 


std of the slope 


% 




Row 


3: 


mean of the std 


% 




Row 


4 : 


std of the std 


% 




Row 


5: 


mean of the skewness 


% 




Row 


6 : 


std of the skewness 


% 




Row 


7: 


mean of the kurtosis 


% 




Row 


8: 


std of the kurtosis 


% 

% 


The 


number 


of columns will be a function of the 



% selected INTERVAL and OVERLAP. 

% 

% If the time series data contains a zero value this will 
% cause errors in the statistics because log(O) = -inf. 

% All beams which contain a zero value are discarded. 

% An error log may be obtained as an output . 
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% Created: May 25, 1999 
% David Schalm 

% 

% Last Modified: May 25, 1999 
% David Schalm 

load data_drive 
drive_num = data_drive ( 1 ) ; 

if nargin == 0 

file_type = input (' Select file type ( sd, cp, cw) . . . ' ) ; 

end 

if file_type == 'sd' 

filepath = { [drive_num . . . 

' : \LWAD 99-1 Data\Tape 17\sd\'] ... 

[drive_num ' : \LWAD 99-1 Data\Tape 18\sd\']}; 
elseif file_type == 'cp' 

filepath = { [drive_num . . . 

' : \LWAD 99-1 DataXTape 17\vd\cp\'] ... 

[drive_num ' : \LWAD 99-1 DataXTape 18\vd\cp\ ' ] } ; 
elseif file_type = = 'cw' 

filepath = { [drive_num . . . 

' : \LWAD 99-1 DataXTape 17\vd\cw\'] ... 

[drive_num ' : \LWAD 99-1 DataXTape 18\vd\cw\ ' ] } ; 

end 

row_index = 0 ; 
error_index = 0 ; 

for b = 1:1:2 
clear fname 

D = dir (char ( filepath (b) )) ; 

[ fname{ 1 : length (D) } ] = deal (D . name) ; 

% get rid of ' . ' and ' . . ' filenames 
fname = fname ( 3 : length ( fname ) ) ; 

for m = 1 : 1 : length ( fname) 

load ( [char (filepath (b) ) char ( fname (m) ) ] ) 
disp ([' Current file is ' filepath{b} fname{m}]) 

X = ( 1 : 1 : length (bmf ))' ; 

N = round (samplerate* interval ) ; 
num_of_beams = size (bmf, 2); 

% Create time vector 
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data_start = 1; 

data_duration = length (bmf ) /samplerate - 
( 1/samplerate) . . . 

+ data_start; 

time = (data_start: 1/samplerate :data_duration) ' ; 
[mtime,ao] = mmean ( time, N, overlap) ; 

for n = 1 : 1 :num_of_beams 

% Skip over data with 0 values --> cause errors 
% in statistics 
if sum(bmf(:,n) == 0) >0 
else 

row_index = row_index + 1; 

% Calculate amplitude of data (dB) 
log_data = 20*logl0 (abs (bmf ( : , n) ) ) ; 

mean_data = mmean (log_data, N, overlap) ; 
lwad_stats_slope (row_index, : ) = ... 
gradient (mean_data) ; 

% Detrend the data 
data = log_data - ... 

polyval (polyf it (X, log_data, 3 ) , X) ; 

lwad_stats_std (row_index, : ) = ... 

ms td (data, N, overlap) ; 
lwad_stats_skewness (row_index, : ) = ... 

mskewness (data, N, overlap) ; 
lwad_stats_kurtosis (row_index, : ) = ... 
mkurtosis (data, 1 , N, overlap) ; 

end 

end 

end 

end 

[M,N] = size (lwad_stats_slope) ; 

lwad_stats_slope ( (M+l) , : ) = mean (lwad_s tat s_s lope) ; 
lwad_s tat s_s lope ( (M+2 ) , : ) = std ( lwad_stats_slope ( 1 :M, : ) ) 

lwad_stats_std ( (M+l) , : ) = mean ( lwad_s tat s_std) ; 
lwad_stats_std ( (M+2 ) , : ) = std(lwad_stats_std(l:M,:)); 

lwad_stats_skewness ( (M+l) , : ) = mean (lwad_s tat s_skewness) 
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lwad_stats_skewness ( (M+2) , : ) = ... 
std ( lwad_stats_skewness (1 :M, :)); 

lwad_stats_kurtosis ( (M+l ) , : ) = mean ( lwad_stats_kurtosis ) 
lwad_stats_kurtosis ( (M+2) , : ) = ... 
std ( lwad_stats_kurtosis ( 1 :M, :)) ; 



area_stats (1:2, : ) 
area_stats (3:4, : ) 
area_stats (5:6, : ) 
area_s tats (7:8, : ) 



lwad_stats_slope ( (M+l) : (M+2) , : ) ; 
lwad_stats_std ( (M+l) : (M+2) , : ) ; 
lwad_stats_skewness ( (M+l) : (M+2) , : ) ; 
lwad_stats_kurtosis ( (M+l) : (M+2) , : ) ; 
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function [y , actual_overlap] = ramean (x, N, overlap) 

%MMEAN Moving mean value . 

% MMEAN(x,N, overlap) should be used to calculate a "moving 
% average" where N is the number of points to use for each 
% interval, and OVERLAP is the ammount of overlap to use. 

% In some cases an assigned overlap may not result in an 

% integer value for STEP (i.e. number of points to advance 

% after each calculation) . In this case, OVERLAP will be 

% rounded to the nearest value which will result in an 

% integer value for STEP. The ACTUAL_OVERLAP can be 
% displayed if desired. 

% 

% See also MEAN, MSTD, MKURTOSIS, MS NEWNESS , LWAD . 

% Created: May 8, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 

if nargin == 2 
overlap = 0 ; 

end 

step = round (( 1-overlap) *N) ; % Number of points to advance 

i = 1; 

while (N+ ( i-1 ) *step) <= length (x) 

y(i) = mean(x( (1+ (i-1) *step) : (N+ (i-1) *step) ) ) ; 
i = i + 1 ; 

end 

if nargout == 2 

actual_overlap = l-(step/N); 

end 
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function [y , actual_overlap] = mstd (x, N, overlap) 

%MSTD Moving standard deviation 
% 

% MSTD (x, N, overlap) should be used to calculate a "moving 
% std" where N is the number of points to use for each 

% interval, and OVERLAP is the ammount of overlap to use. 

% In some cases an assigned overlap may not result in an 

% integer value for STEP (i.e. number of points to advance 

% after each calculation) . In this case, OVERLAP will be 

% rounded to the nearest value which will result in an 

% integer value for STEP. The ACTUAL_OVERLAP can be 
% displayed if desired. 

% 

% See also MEAN, MEDIAN, STD, MIN, MAX, COV, LWAD . 

% Created: May 8, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 

if nargin == 2 
overlap = 0 ; 

end 

step = round ( (1-overlap) *N) ; % Number of points to advance 

i = 1; 

while (N+ ( i- 1 ) *s tep) <= length(x) 

y(i) = std(x( (1+ ( i — 1 ) *step) : (N+ ( i — 1 ) *step) ) ) ; 
i = i + 1 ; 

end 

if nargout == 2 

actual_overlap = l-(step/N); 

end 
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function [ s , actual_overlap] = mskewness (x, N, overlap) 
%MSKEWNESS Moving Skewness . 

% MSKEWNESS (x, N, overlap) should be used to calculate a 
% "moving skewness" where N is the number of points to use 
% for each interval, and OVERLAP is the ammount of overlap 
% to use. In some cases an assigned overlap may not 
result 

% in an integer value for STEP (i.e. number of points to 
% advance after each calculation) . In this case, OVERLAP 

% will be rounded to the nearest value which will result 

% in an integer value for STEP. The ACTUAL_OVERLAP can be 
% displayed if desired. 

% 

% See also SKEWNESS, MMEAN, MSTD, MKURTOSIS, LWAD . 

% Created: May 8, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 

if nargin == 2 
overlap = 0 ; 

end 

step = round (( 1-overlap) *N) ; % Number of points to advance 

i = 1; 

while (N+ ( i-1 ) *step) <= length(x) 

s(i) = skewness (x( (1+ (i-1) *step) : (N+ (i-1) *step) )) ; 
i = i + 1 ; 

end 

if nargout = = 2 

actual_overlap = l-(step/N); 

end 
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function [k, actual_overlap] = mkurtosis (x, flag, N, overlap) 
%MKURTOSIS Moving Kurtosis . 

% MKURTOSIS (x, flag, N, overlap) should be used to calculate 
% a "moving kurtosis" where N is the number of points to 
% use for each interval, and OVERLAP is the ammount of 
% overlap to use. In some cases an assigned overlap may 

% not result in an integer value for STEP (i.e. number of 

% points to advance after each calculation) . In this 
% case, OVERLAP will be rounded to the nearest value 
% which will result in an integer value for STEP. The 
% ACTUAL_OVERLAP can be displayed if desired. 

% 

% MKURTOSIS (x, 1,N, overlap) calculates excess kurtosis. 

% 

% See also KURTOSIS, MSTD, MSKEWNESS, MMEAN, LWAD. 

% Created: May 8, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 

if nargin == 3 
overlap = 0 ; 

end 

step = round (( 1-overlap) *N) ; % Number of points to advance 

i = 1; 

while (N+ (i-1) *step) <= length(x) 
if flag == 1 

k(i) = kurtosis (x( (1+ (i-1) *step) : (N+ (i-1) *step) ) ) - 3; 

else 

k(i) = kurtosis (x( (1+ (i-1) *step) : (N+ (i-1) *step) )) ; 

end 

i = i + 1 ; 

end 

if nargout == 2 

actual_overlap = l-(step/N); 

end 
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function [cglatout, cglonout , tgtlatout , tgtlonout] = ... 

lwad_posit (data_f ile) 

%LWAD_POSIT 

% Determine the coordinates of the DD (lat/lon) based on 
% the ping time . 

% 

% See also LWAD 

% Created: May 15, 1999 

% David Schalm 



% Last Modified: July 26, 1999 
% David Schalm 

load (data_f ile, ' -mat ' ) 
load evt31.mat 

ref_time = s tr2num (ping_time ( 1 : 2 ) ) *3 600 + ... 

str2num(ping_time (4 : 5 ) ) *60 + str2num (ping_time ( 7 : 8 ) ) ; 

i = 1; 

while evt31_time ( i) <= ref_time 
i = i + 1 ; 

end 



interp_factor = (evt31_time(i) - ref_time) /60 ; 

cglatout = cg_lat(i) - (cg_lat(i) - cg_lat ( i-1 ) ) * . . . 
in terp_f ac t or ; 

cglonout = cg_lon(i) - (cg_lon(i) - cg_lon ( i-1 ) ) * . . . 
interp_f actor ; 

tgtlatout = tgt_lat(i) - (tgt_lat(i) - tgt_lat ( i-1 ) ) * 
interp_f actor ; 

tgtlonout = tgt_lon(i) - (tgt_lon(i) - tgt_lon ( i-1 ) ) * 
interp_f actor ; 
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function [course, speed, ctrbrng, botdepth, deangle] = ... 

lwad_prm (data_f ile) 

%LWAD_PRM 

% This function will access the parameter file 
% ' evt31prm.mat ' to determine the ping parameters 

% associated with the selected data file. 



% 

% 

% 

% 

% 



course : 
speed : 
ctrbrng : 
botdepth : 
deangle : 



ships course in degrees True 
ships speed in knots 

center of ping section in degrees True 

bottom depth in fathoms 

sonar depression angle in degrees 



% See also LWAD 



% Created: May 11, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 



load (data_f ile, ' -mat ' ) 
load evt31prm.mat 

ref_time = str2num (ping_time (1 : 2 ) ) *3600 + ... 

str2num(ping_time (4 : 5 ) ) *60 + str2num (ping_time (7 : 8 ) ) ; 

i = 1; 

while timeday(i) < ref_time 
i = i + 1 ; 

end 

i = i-1 ; 

course = coursed); 
speed = speed ( i ) ; 
ctrbrng = ctrbrng (i); 
botdepth = bo tdepth (i); 
deangle = deangle (i); 
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function lwad_map ( fname) 
%LWAD MAP 



% 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



% 

% 

% 



Create a map of the exercise area 
(1) A map of the LWAD 99-1 exercise will be created 
showing the following: 

(a) blue eg track history, red circle is eg 
current location 

(b) red tgt track history, blue circle is tgt 
current location 

(c) ping section is shown in yellow with dashed 
yellow center section true bearing line. 

Yellow asterisks show lOnm range. 

(d) latitude/ longitude grid lines (black dashed 
lines) 

(e) Sediment map showing hal hook (black circles) 
with phi size approx 4-7, clay silt (black dots) 
with phi size approx 6-8, and sand 

(black slashed lines) with phi size approx 1. 
This map is based on data received from 
J. Fulford (NRL Stennis) . 



If no filename is specified, 
the user to select a file. 



a gui will open to allow 



See also LWAD 



% Created: May 15, 1999 
% David Schalm 

% 

% Last Modified: July 26, 1999 
% David Schalm 

if nargin == 0 

load data_drive 

[ filename , pathname] = uigetf ile ( [data_drive] ) ; 
fname = [pathname filename] ; 
data_drive = pathname; 

save 'e:\schalm\lwad 99-l\data_drive .mat ' data_drive 

end 



load ( fname , ' -mat ' ) 

[eglat, eglon, tgtlat, tgtlon] = lwad_posit ( fname) ; 

[course, speed, ctrTbrng, bo tdepth, deangle] = lwad_prm ( fname ) ; 
load evt31.mat; 
num_of_beams = size(bmf,2); 

figure (100) 
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elf 

set ( 100 , . . . 

'NumberTitle ' , 'off', ... 

' Name ' , ' LWAD 99-1 Event 31'); 

% Plot the sediment 

pic = imread (' f lor ida' , 'bmp' ) ; 

sed_lon = linspace ( -84 . 5 , -83 . 9 , 600 ) ; 

sed_lat = linspace (25 . 4 , 25 . 9 , 600) ; 

image (sed_lon, sed_lat, pic) ; 

hold on, grid on,’ axis xy 

%plot eg and tgt track history- 
plot ( -cg_lon, cg_lat) , hold on 

xlabel (' Longitude (W) ' ) , ylabel (' Latitude (N) ' ) 

title ( [ ' LWAD 99-1 Event 31']), axis equal, axis image 

plot ( -eglon, eglat , 'ro' ) 

plot ( -tgt_lon, tgt_lat , 'r') 

plot (-tgtlon, tgtlat, 'bo' ) 

% Plot a lOnm range circle, and ctr bearing line 
theta_one = ctrTbrng - num_of_beams*5/2 ; 
theta = theta_one:5:(theta_one+(num_of_beams-l)*5); 
theta (num_of_beams+l ) = theta_one; 

for i = 1 : 1 : (num_of_beams+l ) 

circle_lat ( i) = eglat + 10*cos (theta(i) *pi/180) /60 
circle_lon ( i ) = -(eglon - 10*sin(theta(i) *pi/180) / 
(60*cos (cglat*pi/180 ) ) ) ; 

end 

if num_of_beams == 24 

plot (circle_lon, circle_lat , 'y* ' ) 
plot ( [-eglon circle_lon ( 13 ) ] , . . . 

[eglat circle_lat (13) ] , ' --y' ) 
circle_lat = [eglat circle_lat ( 1 : 24 ) eglat]; 
circle_lon = [-eglon circle_lon ( 1 : 24 ) -eglon]; 
plot (circle_lon, circle_lat, 'y' ) 
else 

plot (circle_lon, circle_lat, 'y* ' ) 

plot ([-eglon circle_lon ( (num_of_beams/2 ) + 1)],... 

[eglat circle_lat ( (num_of_beams/2 ) + 1)],' — y' ) 
plot (circle_lon, circle_lat , 'y' ) 
end 
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function [cout] = lwa d_con tour ( f name , interval , overlap) 

% LWAD_C ONTOUR 

% Create a sediment contour based on the selected ping. 

% 

% This plot is created by determining the difference 
% between area averaged slope data and slope of 

% individual beam data. 

% 

% If no filename, interval, or overlap are specified, a 

% gui will open to allow the user to select a file. 

% The interval and overlap will default to 1 sec and 

% .75. Modify this file in order to change the 

% defaults . 

% 

% The plot shows the eg current position (black circle) , 
% target current position (red circle) , and exercise 

% area (black dashed lines) . 

% 

% This function uses code from the Matlab function 

% ' contour. m' 

% 

% See also LWAD_CONTOURF , CONTOUR, LWAD. 

Created: May 15, 1999 
David Schalm 

Last Modified: July 26, 1999 
David Schalm 



if nargin == 0 

load data_drive 

[ filename , pathname] = uigetf ile ( [data_drive '*.*']); 
fname = [pathname filename] ; 
data_drive = pathname; 

save 'e:\schalm\lwad 99-l\data_drive.mat' data_drive 

%%%%%% Default values %%%%% 
interval = 1 ; 
overlap = .75; 

spacing = interval - interval* overlap; 

%interval = input ('What is the interval (.5,1,2): '); 

%overlap = input ('What is the overlap (0,. 25,. 5,. 75): '); 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 



else 

if fname ( (length (fname) -2 ): length ( fname) ) == 'mat' 

filename = fname ( (length (fname) -23 ): length ( fname) ) ; 
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else 

filename = fname (( length ( fname) -19 ): length ( fname) ) ; 

end 

end 

% Load the data file 
load ( fname , ' -mat ' ) 

[cglat , cglon, tgtlat , tgtlon] = lwad_posit ( fname) ; 

[ course , speed , c t rTbrng , bo tdep th , deangle ] = lwad_prm ( fname ) ; 

if f ilename ( 5 : 6 ) == 'SD' | f ilename ( 5 : 6 ) == 'sd' 
file_type = filename (5 : 6 ) ; 
else 

file_type = filename (19 : 20) ; 

end 

% Load the average statistics 

area_stats_fname = ['e:\schalm\lwad 99-l\area_stats\ ' ... 

file_type '_stats' num2str ( interval) '.mat']; 
load (area_stats_f name) 
if overlap == 0 

area_stats = stats_0; 
elseif overlap = = .25 

area_stats = stats_25; 
elseif overlap == .5 

area_stats = stats_5; 
else 

area_stats = stats_75; 

end 

% Create time\range vector 
sound_speed = .81; %nm/s 
N = round ( interval* samplerate ) ; 
data_start = 1 ; 

data_duration = length (bmf) /samplerate - ( 1/ samplerate) .. . 

+ data_start; 

time = (data_start : 1/ samplerate : data_durat ion) ' ; 

[mtime,ao] = mmean ( time, N, overlap) ; 
mrange = sound_speed*mtime/2 ; 

% Create bearing vector 
num_of_beams = size(bmf,2); 

theta_one = ctrTbrng - num_of_beams*5/2 + 5/2; 
theta = theta_one : 5 : ( theta_one+ (num_of_beams-l) *5 ) ; 

index = 0 ; 

for i = 1 : 1 : num_of_beams 
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if sum (bmf ( : , i ) == 0) >0 
% skip this beam 
else 

log_data = 20*logl0 (abs (bmf ( : , i) ) ) ; 

index = index + 1; 

thetal ( index) = theta (i) ; 

mean_data = mmean ( log_data , N , overlap ) ; 

slope_data = (gradient (mean_data) ) ; 

dif f_data ( index, : ) = slope_data - area_stats (1 , : ) ; 

end 

end 

cl = contourc (mrange, thetal , dif f_data) ; 

i = 1; 

while i < length (cl) 
nexti = i + cl(2,i); 
c_new ( : , i ) = c 1 ( : , i ) ; 

% convert range/brng to lat/lon 
for k = ( i+1) : 1 : nexti 

c_new(l,k) = -(cglon - sin (cl (2 , k) *pi/180 ) * . . . 

cl ( 1 , k) / ( 60*cos (cglat*pi/180 ) ) ) ; 
c_new(2,k) = cglat + cos (cl (2, k) *pi/180) * . . . 
cl (l,k) /60; 

end 

i = nexti + 1; 

end 

c = c_new; 
figure 

lims = [-84.4 -84 25.5 25.8]; 

colortab = get (newplot, ' colororder ' ) ; 

[mc,nc] = size (colortab) ; 

% Create the plot 
newplot ; 
if -ishold 

view(3); grid on 

set (gca, 'xlim' , lims (1:2) , 'ylim' , lims (3:4) ) 
end 

limit = size(c,2); 

i = 1; 
h = [ ] ; 
color_h = [ ] ; 
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while(i < limit) 

z_level = c(l,i); 
npoints = c(2,i); 
nexti = i+npoints+1; 

xdata = c (1, i+1: i+npoints) ; 
ydata = c (2 , i+1 : i+npoints) ; 

zdata = z_level + 0*xdata; % Make zdata the same size 

% as xdata 

% Create the patches or lines 

cu = patch ( 'XData' , [xdata NaN] , 'YData' , [ydata NaN] , . . . 

'ZData' , [zdata NaN] , 'CData' , [zdata NaN] , ... 

' f acecolor ' , ' none ' , ' edgecolor ' , ' flat ' , . . . 

'userdata' , z_level) ; 
h = [h; cu(:) ] ; 

color_h = [color_h ; z_level] ; 
i = nexti; 

end 

for i = 1: length (h) 

set (h (i) ; ' Zdata' , [] ) ; 
end 

view (2 ) ; 

set (gca, 'box' , 'on') ; 

axis equal, grid on, zoom on, hold on 
axis ( [-84.5 -83.9 25.4 25.9]) 

xlabel (' Longitude (VI)'), ylabel (' Latitude (N) ' ) 
warning off 

title ([' Event 31/' f ilename ( 13 : 16 ) ', ' ... 

ping_date ' , ' ping_time] ) 

caxis ( [-12 12] ) 
colorbar 

plot ( -cglon, cglat , 'ko' ) 
plot (-tgtlon, tgtlat , 'ro') 

% Plot eg box 

plot ( [-84 . 27 -84.27 -84.11 -84.11 -84.27],... 

[25.61 25.71 25.71 25.61 25.61], 'k—') 
zoom on 
brighten (1) 

if nargout > 0 
cout = c ; 

end 



94 



oV oV c 69 dP 



function [cout] = lwad_contourf ( f name , interval , overlap) 
%LWAD_CONTOURF 

% Create a filled sediment contour based on the selected 
% ping . 

% 

% This plot is created by determining the difference 
% between area averaged slope data and slope of 

% individual beam data. 

% 

% If no filename, interval, or overlap are specified, a 
% gui will open to allow the user to select a file. 

% The interval and overlap will default to 1 sec and 

% .75. Modify this file in order to change the 

% defaults. 

% 

% The plot shows the eg current position (black circle) , 

% target current position (red circle) , and exercise 

% area (black dashed lines) . 

% 

% This function uses code from the Matlab function 
% 'contourf.m' 

% 

% See also LWAD_CONTOUR, CONTOURF, LWAD. 

Created: May 15, 1999 
David Schalm 

Last Modified: July 26, 1999 
David Schalm 

if nargin == 0 

load data_drive 

[filename, pathname] = uigetf ile ( [data_drive '*.*']); 
fname = [pathname filename] ; 
data_drive = pathname; 

save 'e:\schalm\lwad 99-l\data_drive .mat ' data_drive 

%%%%%% Default values %%%%% 
interval = 1; 
overlap = .75; 

%interval = input ('What is the interval (.5,1,2): '); 

%overlap = input ('What is the overlap (0,. 25,. 5,. 75): ') 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

else 

if fname (( length ( fname) -2 ): length ( fname) ) == 'mat' 

filename = fname (( length ( fname) -23 ): length (fname)); 
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else 

filename = fname ( (length (fname) -19) : length (fname) ) ; 

end 

end 

% Load the data file 
load ( fname , ' -mat ' ) 

[cglat , cglon, tgtlat , tgtlon] = lwad_posit ( fname) ; 

[course, speed, ctrTbrng,botdepth, deangle] = lwad_prm ( fname) ; 
spacing = interval - interval* overlap; 

if filename (5 : 6) == 'SD' | filename (5 : 6 ) == 'sd' 
file_type = filename (5 : 6 ) ; 
else 

file_type = filename (19 :20) ; 

end 

% Load the average statistics 

area_s tat s_f name = ['e:\schalm\lwad 99-l\area_stats\ ' ... 

file_type '_stats' num2str (interval) '.mat']; 
load ( area_s tats_f name ) 
if overlap == 0 

area_stats = stats_0; 
elseif overlap == .25 

area_stats = stats_25; 
elseif overlap == .5 

area_stats = stats_5; 
else 

area_stats = stats_75; 

end 

% Create time\range vector 
sound_speed = .81; %nm/s 
N = round ( interval *samplerate) ; 
data_start = 1; 

data_duration = length (bmf ) /samplerate - (1/samplerate) . . . 

+ data_start; 

time = (data_start : 1/samplerate : data_duration) ' ; 

[mtime,ao] = mmean ( t ime , N , overlap) ; 
mrange = sound_speed*mtime/2 ; 

% Create bearing vector 
num_of_beams = size (bmf, 2); 

theta_one = ctrTbrng - num_of_beams*5/2 + 5/2; 
theta = theta_one: 5 : (theta_one+ (num_of_beams-l) *5) ; 
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index = 0 ; 

for i = 1 : 1 : num_of_beams 

if sum (bmf ( : , i) == 0) >0 
% skip beams w/ bad data 
else 

log_data = 20*logl0 (abs (bmf ( : , i) ) ) ; 

index = index + 1 ; 

thetal (index) = theta (i); 

mean_data = mmean ( log_data, N, overlap) ; 

slope_data = (gradient (mean_data) ) ; 

diff_data ( index, : ) = (slope_data - area_stats ( 1 , : ) ) ; 

end 

end 

diff_data = dif f_data/ spacing ; 

lin = ' ' ; 
col = ' ' ; 
x = mrange; 
y = thetal; 
z = diff_data; 
nv = [ ] ; 

%nv = 4 ; 

if (size (y, 1) ==1) , y=y' ; end; 
if (size (x, 2) ==1) , x=x'; end; 

[mz , nz] = size ( z ) ; 

lims = [-84.4 -84 25.5 25.8]; 

i = find(isfinite(z) ) ; 
minz = min ( z ( i ) ) ; 
maxz = max ( z ( i ) ) ; 

if length (nv) <= 1 
if isempty(nv) 

CS=contourc ( [minz maxz ; minz maxz] ) ; 
else 

CS=contourc ( [minz maxz ; minz maxz],nv); 
end 

% Find the levels 
ii = 1; 

nv = minz; % Include minz so that the contours are 
%totally filled 
while (ii < size(CS,2)), 
nv= [nv CS ( 1 , ii ) ] ; 
ii = ii + CS(2,ii) + 1; 
end 
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end 



% Don't fill contours below the lowest level specified in 
% nv. To fill all contours, specify a value of nv lower 
% than the minimum of the surface. 
draw_min=0 ; 
if any(nv <= minz) , 
draw_min=l ; 
end 

% Get the unique levels 
nv = sort ( [minz nv (:)']); 
zi = [1, f ind(dif f (nv) ) +1] ; 
n v = n v ( z i ) ; 

% Surround the matrix by a very low region to get closed 
% contours, and replace any NaN with low numbers as well. 

zz=[ repmat (NaN, 1 , nz+2 ) ; repmat (NaN,mz, 1) z ... 

repmat (NaN, mz , 1 ) ; repmat (NaN, 1 , nz+2 )] ; 

kk=f ind ( isnan ( zz ( : ) ) ) ; 

zz (kk) =minz-le4* (maxz-minz) +zeros (size (kk) ) ; 

xx = [2*x ( : , 1) -x ( : , 2 ) , x, 2 *x ( : , nz ) -x ( : , nz-1 ) ] ; 
yy = [2*y (1, : ) -y (2, : ) ; y; 2 *y (mz , : ) -y (mz-1 , : ) ] ; 
if (min (size (yy) ) ==1) , 

[CS,msg] =contours (xx, yy, zz , nv) ; 
else 

[CS,msg] =contours (xx ( [ 1 l:mz mz] , : ) ,yy ( : , [1 l:nz nz] ) . . . 

, z z , nv ) ; 

end; 

if -isempty (msg) , error (msg); end 
c_new = [ ] ; 

i = 1; 

while i < length (CS) 

nexti = i + CS(2,i); 
c_new(:,i) = CS(:,i); 

% convert range/brng to lat/lon 
for k = (i+1) : l:nexti 

c_new(l,k) = - (cglon - sin (CS (2 , k) *pi/180 ) * . . . 

CS (l,k) / (60*cos (cglat*pi/180 ) ) ) ; 
c_new(2,k) = cglat + cos (CS (2 , k) *pi/180 ) * . . . 

CS (l,k) / 60 ; 

end 

i = nexti + 1; 
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end 

CS = c_new; 
figure 

% Find the indices of the curves in the c matrix, and get 
% the area of closed curves in order to draw patches 
% correctly, 
ii = 1; 
ncurves = 0 ; 

I = []; 

Area= [ ] ; 

while (ii < size(CS,2)), 
nl=CS (2 , ii ) ; 
ncurves = ncurves + 1; 

I (ncurves) = ii; 

xp=CS ( 1 , ii+ ( 1 : nl ) ) ; % First patch 

yp=CS (2 , ii+ (1 : nl ) ) ; 

Area (ncurves ) =sum ( dif f (xp) . * (yp ( 1 : nl-1 ) +yp ( 2 : nl ) ) /2 ); 
ii = ii + nl + 1; 
end 

newplot; 
if ~ishold, 
view (2) ; 

set (gca, 'xlim' , lims (1:2) , 'ylim' , lims (3:4) ) 
end 

% Plot patches in order of decreasing size. This makes sure 
% that all the leves get drawn, not matter if we are going 
% up a hill or down into a hole. When going down we shift 
% levels though, you can tell whether we are going up or 
% down by checking the sign of the area (since curves are 
% oriented so that the high side is always the same side) . 

% Lowest curve is largest and encloses higher data always. 

H= [ ] ; 

[FA, IA] =sort ( -abs (Area) ) ; 
if -isstr (get (gca ,' color ')) , 
bg = get (gca ,' color ') ; 
else 

bg = get (gcf ,' color ') ; 
end 

if isempty(col) 

edgec = get (gcf , ' defaultsurf aceedgecolor ' ) ; 
else 

edgec = col; 
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end 

if isempty(lin) 

edgestyle = get (gcf , ' defaultpatchlinestyle' ) ; 
else 

edgestyle = lin; 
end 

for j j =IA, 

nl=CS (2, I ( j j ) ) ; 
lev=CS (1, 1 ( j j ) ) ; 
if (lev ~=minz | draw_min ) , 
xp=CS ( 1 , 1 ( j j ) + (1 :nl) ) ; 
yp=CS (2 , I ( j j ) + (1 :nl) ) ; 

if (sign(Area(j j ) ) ~=sign (Area (IA(1) ) ) ), 

kk=f ind (nv==lev) ; 

if (kk>l+sum (nv<=minz ) * (~draw_min) ) , 
lev=nv(kk-l) ; 
else 

lev=NaN; % missing data section 

end; 
end; 

if (isf inite (lev) ) , 

H= [H; patch (xp,yp, lev, ' facecolor' , ' flat' , . . . 

' edgecolor ' , edgec , ' linestyle ' , edgestyle) ] 

else 

H= [H ; patch (xp, yp, lev, ' facecolor' , bg, ... 

' edgecolor' , edgec, ' linestyle' , edgestyle) ] 

end; 

end; 

end; 

numPatches = length (H); 
if numPatches>l 

for i=l : numPatches 

set(H(i), ' faceoffsetf actor ' , 0, ... 

' f aceof f setbias ' , ( le-3 ) + (numPatches-i) / . . . 

(numPatches- 1) /30) ; 

end 

end 

% Set the plot options 

axis equal, grid on, zoom on, hold on 

axis ( [-84.5 -83.9 25.4 25.9]) 

xlabel (' Longitude (W) ' ) , ylabel (' Latitude (N) ' ) 
warning off 

title ([' Ping : ' f ilename ( 13 : 16 ) ... 

',Mode: SD, Interval: 1, Overlap: .75']) 
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% Uncomment the next line for a b/w plot 
%colormap (gray) 

caxis ( [-12 12]) 
colorbar 

plot (-cglon, cglat, 'ko' ) 
plot (-tgt Ion, tgtlat , ' ro' ) 

% Plot eg box 

plot ( [-84.27 -84.27 -84.11 -84.11 -84.27],. 

[25.61 25.71 25.71 25.61 25.61], 'k—') 
zoom on 

if nargout > 0 
cout = c ; 

end 
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% COMPOSIT 
% 

% This script file will create a composite graph of many 
% separate pings. For now, it will only work for SD mode 
% pings. 

% Location of the pings to combine 

filepath = {'e:\schalm\LWAD 99-1 Data\Tape 17\sd\'}; 

%%%% Default values (Modify if desired) %%%% 
interval = 1 ; 
overlap = .75; 

spacing = interval - interval*overlap; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Load the area stats file 

area_stats_fname = ['e:\schalm\lwad 99-l\area_stats\ ' . . . 
'sd_stats' num2str ( interval) '.mat']; 

load ( area_s tats_f name ) 

if overlap == 0 

area_stats = stats_0; 
elseif overlap == .25 

area_stats = stats_25; 
elseif overlap == . 5 

area_stats = stats_5; 
elseif overlap == .75 

area_stats = stats_75; 

end 

% Get a directory of the pings to combine 
D = dir (char ( filepath) ) ; 

[fname{l : length (D) } ] = deal (D. name) ; 

% get rid of ' . ' and ' . . ' filenames 
fname = f name (3 ; length ( f name) ) ; 

% Select the Lat/Lon you want to look at 
% and the grid resolution 
Ion = -84.45: .01:-83 .9; 
lat = 25.45: .01:25.9; 

compos it_graph = zeros ( length ( Ion) , length ( lat )) ; 
ref_matrix = zeros ( length ( Ion) , length ( lat )) ; 

% Fill in the values of compos it_graph one ping at a time 
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for m = 1:1:1 

%for m = 1 : 1 : length ( f name) 

load( [char ( f ilepath) char ( fname (m) ) ] , ' -mat' ) 

disp ([' Current file is ' char ( f ilepath) char ( fname {m} )] ) 

[cglat, cglon, tgtlat , tgtlon] = lwad_posit . . . 

( [char ( f ilepath) char ( fname (m))]); 

[course, speed, ctrTbrng, botdepth, deangle] = ... 

lwad_prm( [char ( f ilepath) char ( fname (m) ) ] ) ; 
cglon = -cglon; 

for lon_index = 1 : 1 : length ( Ion) 
for lat_index = 1 : 1 : length (lat) 

% check to see if lat/lon is in range of the eg 
dx = 60*cos (cglat*pi/180) * ( Ion ( lon_index) -cglon) ; 
dy = 60* (lat (lat_index) - cglat); 
range = sqrt (dx^2 + dy^2); 

if range > 9.9 

% do nothing because this is out of range 
% of the eg 
else 

brng = atan2 (dy , dx) ; 

% Convert bearings to true (north up) bearings 
if brng < 0 

brngl = -1* (brng*180/pi-90 ) ; 
else 

brngl = 90-brng*180/pi ; 

end 

if brngl < 0 

brngl = brngl+360; 

end 

brngl; 

beam_num = ceil (brngl/5 ) ; 

% Skip over data with 0 values --> cause errors 
% in statistics 

if sum(bmf ( : ,beam_num) == 0) >0 
%do nothing 

% Skip over beams that are +/- 45 degrees 
% of the stern 

elseif course > 075 & course < 105 & ... 
beam_num > 45 & beam_num < 63 
% do nothing 

elseif course > 165 & course < 195 & ... 
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( beam_num < 9 | beam_num > 63) 

% do nothing 
else 

N = round ( interval *samplerate) ; 
log_data = 20*logl0 (abs (bmf ( : ,beain_nuin) ) ) 
mean_data = mmean ( log_data, N, overlap) ; 
slope_data = (gradient (mean_data) ) ; 
diff_data = (slope_data - area_stats . . . 

(1/ : ) ) /spacing; 
sound_speed = .81; %nm/s 
data_time = range*2/sound_speed; 
index = round ( (data_time- . 25 ) / . 25 ) ; 
if index < 1 
index = 1; 

end 

diff_point = diff_data (index) ; 
compos it_graph ( lon_index, lat_index) = ... 
compos it_graph (lon_index, lat_index) + . . 
dif f_point ; 

ref_matrix ( lon_index, lat_index) = ... 
ref_matrix ( lon_index, lat_index) + 1; 

end 

end 

end 

end 

end 

% Divide the lat/lon matrix values by the number of 
% entries at each point 

compos it graph = compos it_graph. /ref_matrix; 

% Plot the composite graph 
figure 

contour f (Ion, la t, compos it_graph' ) , hold on 

% Uncomment the next line for a b/w plot 
%colormap (gray) 
caxis ( [ -8 8 ] ) , 
grid on,colorbar 

xlabel ( 'Longitude (W) ' ) , ylabel ( 'Latitude (N) ' ) 
title ([' Composite Plot']) 

% Plot eg box 

plot ( [-84.27 -84.27 -84.11 -84.11 -84.27],... 

[25.61 25.71 25.71 25.61 25.61], 'k— ') 
zoom on 
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%LWAD_FILE 

% 

% Use this script to open and load an LWAD 99-1 file 
load data_drive 

[filename pathname] = uigetf ile ( [data_drive] ) ; 
fname = [pathname filename] ; 
data_drive = pathname; 

save 'e:\schalm\lwad 99-l\data_drive .mat ' data_drive 

% Load the data file 
load ( fname , ' -mat ' ) 
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